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ABSTRACT 


A program has been written for the identification of parameters 
in certain linear systems. These systems appear in biomedical prob- 
lems, particularly in compartmental models of pharmacokinetics. The 
method presented here assumes that some of the state variables are 
regularly modified by jump conditions. This simulates administration 
of drugs following some prescribed drug regime. 

Parameters are identified by a least-square fit of the linear 
differential system to a set of experimental observations. The method 
is especially suited when the interval of observation of the system is 
very long. 



A PROGRAM FOR IDENTIFICATION OF 
LINEAR SYSTEMS 


1. Computational Methods 

A program for identification of .linear systems is presented. 

By system identification is meant the evaluation of some of the internal 
parameters of a system so that a set of experimental observations on 
the state variables is optimally approximated. 

The type of model considered here is suggested by studies of 
the time history of a drug in the human system using linear compart- 
mental models [ 1] . In this particular field, system identification is an 
important tool not only for the study of the general pharmacokinetic 
behavior of a drug but also for the determination of individual differ- 
ences between subjects. This may help in the prescription of a drug 
regime especially suited to individual characteristics. 

The method used here, expecially useful when the period of 
observation is long, was derived and described in a previous paper [2]. 
The linearity of the system is used to obtain a series of recurrence 
relations that help reduce the amount of computation. 

It is assumed that doses are administered at equal intervals 
and that every dose (including the initial dose) instantaneously increases 
the concentration of the compartment into which it is injected. Between 
doses, it is assumed that the system is described by a linear, homo- 
geneous, vector-matrix differential equation [3]: 

(1) x = Ax . 

The matrix A depends on various parameters a* , ..., a. m 
whose precise values are not known. It is desired to estimate those 
values on the basis of dynamic measurements made on the system just 
before each dose is administered. 
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It is assumed here that the time scale is chosen so that the 
interval A between dosages is 

(2) A = 1 

The injection of a dose at times t = 1, . . . , N - 1 is simulated 
by the jump condition 

(3) x(i + 0) - x(i - 0) = d , i = 1, 2, . . . , N - 1 

The initial condition is 

(4) x(0) = c . 

At times t 4 i (i = 1, 2, . . . , N) , the state vector x is controlled by 
the relation of Equation (1). The first k variables of x; x\ , x^ , 
are assumed to have been observed at times t = i, i = 1, . . . , N . Let 
the values of these observations be 


(5) 


b J (i) 


j = 1 , . . . , k; i = 1 , 


N , 


where b^(i) stands for the observation of the j-th variable at the time 
t = i . 

The method requires knowledge of an approximation , . . . , 

a^° ) for the desired parameters. Proceeding from this approximation 
and using a combination of a Newton-Raphson scheme and a series of 
recurrence relations, an iterative scheme is followed which eventually 
minimizes the quadratic form 


S = 



Y. (x J (i) - b j (i) ) 2 


( 6 ) 


> 
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where the Y are appropriate weights, provided that a sufficiently close 
initial approximation is given. 

The method proceeds according to the following steps: 

(i) The following initial value matrix-matrix differential systems 


(7) 


~ 0 = A0 , 0St<l , 

(8) 


0(0) = i , 

where 

I is the unit matrix of dimension r (r = dim x ); 

(9) 


4~ f. = AT. + A.0 , 0 ^ t £ 1 , 

dt J J J 

(10) 


T.(0) = 0 , j = 1, m , 

where 

j 

is an r X r matrix and 

(11) 


A. = S/3a.A , j = 1 , • • . , m ; 
J J 

and 



(12) 


4~ A . =AA. +AT.+AT +A.0 , O^t^l 

dt jq jq q j j q jq 

(13) 


A. (0) = 0 , j,q = 1, , m ; 

where 

ft. 

jq 

are r X r matrices and 

(14) 


2 . 

A, = 5 /dx.BxA , j,q=l,...,m ; 


are integrated from t = 0 to t = 1 . Notice that 
(15) (t) ^ A .(t) , 0 ^ t * 1 . 

jq qj 
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The program integrates the above system of — r (m + 1) (m + 2) 

6 

differential equations by a combination of a Runge-Kutta and Adams - 

2 

Moulton schemes. During the integration a "total of 5r (m-bl)(rri + 2) 

places are required to store auxiliary quantities. When the values of the 

matrices <t> ( 1), y.(l), A, (1), j,q=l,...,m are obtained* the program 

J 1 2 

requires their storage which occupies (m + l)(m + 2) memory locations. 

No further auxiliary values are stored. In the previous description it is 
assumed that A, A., and A . are computed using the current approxi- 

j qj 

mation of a . 

(i 1 ) Denoting by (Ct^ , . .., cc ^ ) the current approximation of 

1 m 

the parameters to be identified, the following quantities are computed 


^ - bia,) ^ 
i=l j=l n 


, n = 1 , . . . , m , 


N k r . ,2 j.. 

2 2 y 3 r <xJ(i) - b,,i,) ^ 


i=l j=l 


d x^ (i ) 


dx J (i) 
9 a 


requiring a total of m(m + 3)/2 memory locations. 
The values of the vectors x(i) , 


x (i ) = d x (i ) / d a , n=l,...,m , 

n n 


x (i) = d x(i)/ da da , n, q = 1 , 
nq n q 


. , m 
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are computed recursively using the relations 

(20) x(i + 1) = 0 (1) x(i + 0) , 

(21) x(i + 0) = x(i) + d , i = 1, . . . , N - 1 , 

(22) x(0) = c , 

(23) x (i + l) = 0(l)x (i)+¥ ( 1 ) x (i + 0 ) , 

n n n 

(24) x (0) =0 , n = 1, . . . , m; i = 1, N - 1 , 

n 

(25) x (i + l) = 0(l)x (i) + Y (1 )x (i) + Y (1 )x (i) + A (l)x(i + 0) 

nq nq q n n q nq 

(26) x (0) = 0 , n, q = 1, . . . , m; i = 1 , . . . , N - 1 . 

nq 

Values of the quantities produced by the recurrent process need 
not be sequentially stored and the quantities defined by Equations (16) 
and (17) may be computed as running sums. The total storage require- 
ment for this step is then — r (m + l)(m + 2) + ^-m (m + 3) . 

(iii) The linear algebraic system 


(27) 


E + 

n 


2 


B (a 1 - a° ) = 0 , n 
nq q q 


1 , . . . , m , 


is solved for & , q = 1, • . • , m . The result is the new approximation 

of the parameter vector a . 
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Steps (i), (ii), and (iii) are repeated until 
(28) | a 1 - a° ! < e , q = 1, . . . , m , 

q q * 

where e is user provided. 

2. Flowchart 

The flowchart in Figure 1 is explained by the description of the 
computational method in the previous section. 

3. Typical Sample Run 

The output of the sample run below includes the following infor- 
mation after the run number identification. 

(a) The input constants IC(i), i = 1, . . . , 13 which are user 
supplied. They are integer constants of the program and integer control 
parameters. Their function is described in the software specification 
section. 

(b) The input constants RC(i) which are noninteger constants 
to be supplied by the user. 

(c) The initial conditions x(0) = c . 

(d) The initial approximations for the parameters to be iden- 
0 

tified (a , q = 1, . . . , m ). 

q 

(e) The drug dosages used during the observation period. 

(f) The experimental observations on the system. 

From this point on, the following output is printed at each step of the 
iterative process. 

(g) The iteration cycle number. 

(h) The current value of the parameters to be identified. 
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(i) Optional. The value of the matrix 0 (t) , as computed by 
the integration, is printed every ]V< P steps* This option was used in 
the sample run with MP = 100 . 

(j) Optional. The values x (i - 0) are printed at every obser- 
vation point* This option was not used in the sample run. 

(k) Optional. The values B are printed in matrix form B 
and the values E^ in vector form E . This option was used in the 
sample run. 

(l) Optional. Prints the matrix B * , the value of the deter- 
minant of B \ the vector B ^E , and the matrix BB * This option 
is used to check the reliability of the matrix inversion and was used in 
the sample run. 

(m) The corrected value of the parameters* 

When the iteration ends, a message is printed letting the user know 
whether convergence has occurred or the allowed number of iterations 
has been exceeded. 

The values of x(j - 0) at the observation points using the identi- 
fication parameters are printed. 

The particular identification problem studied in the sample run 
is the determination of the constants u , u in the system. 

1 u 

*1 (t) = x 2 (t) 

x 2 (t) = - (t) + u 2 x 2 (t) , i s: t < i + 1 , i = 0, 1, ...» 9 . 

x^ (i +0) - x^ (i - 0) = 0 
x (i + 0) - x 2 (i - 0) = 1 , i = 1 , . . . , 10 

with the initial conditions 
X L (0) = 0 


V°) = 1 
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Observations of the variable x^ used in the sample run were 
generated at times t = 1, 2, .,.,10 by exact integration of the above 
system using the values u^ = 10 and u^ = 11 . No observations of the 
variable x^ were made available to the program. 

In this example, the solutions can be expressed in the form 

X 1 (t ) = A^” t + A^e " 10t 
x 2 (t) = B^" t + B 2 e " 10t 

where the numbers A A , B , B depend on the particular value 

J ■ Cj i. L* 

(i, i + 1) to which t belongs but are otherwise independent of time. 

Since e * is much larger than e except in a neighborhood 

of zero, this seemingly simple differential system creates a difficult 
identification problem. Initial estimates for the parameters u and u 

X Lr 

were = 8 , = 8 . 
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Sample Run 


INPUT CONSTANTS IC(I) 


RUN NUMBER 


1 

2 

2 

2 

3 

2 

4 

1 

5 

10 

6 

10 

7 

0 

8 

1 

9 

0 

10 

1 

11 

1 

12 

100 

13 

3 

INPUT CONSTANTS RC(I) 

1 

0. 99999979E-02 

2 

0. 99999993E-03 

3 

0. 10000000E-01 

INITIAL CONDITIONS C(I) 

1 

0 . 0 

2 

0. 1 0000000E-01 

INITIAL APPROXIMATION 

1 

0. 80000000E-01 

2 

0. 80000000E-01 

DOSAGES 

1.0 

0. 100000E-01 

2.-.0 

0. 100000E-01 

3. 0 

0. 1 00000E-01 

4. 0 

0. 100000E-01 

5.0 

0. 100000E-01 

6.0 

0. 100000E-01 

7.0 

0. 100000E-01 

8. 0 

0. 100000E-01 

9. 0 

0. 100000E-01 

10. 0 

0. 100000E-01 
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OBSERVATIONS 

0. 40870443E-01 
0. 55907674E-01 
0. 61439544E-01 
0. 63474536E-01 
0. 64223170E-01 
0. 64498603E-01 
0. 64599931E-01 
0. 64637244E-01 
0. 64650953E-01 
0. 64655960E-01 

CYC LE NUMBER 1 

OLD CONSTANTS 

UC( 1)= 0. 80000000E-01 

UC( 2) = 0. 80000000E-01 

INTEGRATION OF PHI AT TIME = 1. 000 
0. 373832E-00 0.545867E-01 

-0. 436703E-00 -0. 628696E-01 

SECOND PARTIALS OF S MATRIX 

0.244998E-02 -0. 4101 95E-03 

-0.410195E-03 0. 994315E-05 

PARTIAL OF S VECTOR 

-0. 203786E-02 
0. 444756E-03 

INVERSE OF SECOND PARTIALS OF S MATRIX 

-0. 690977E-02 -0. 285057E -04 

-0. 285057E-04 -0. 170256E-05 

DETERMINATE = -0. 143900E-06 

PRODUCT OF MATRIX INVERSE AND PARTIAL VECTOR 

-0. 112699E-01 
-0. 176319E-01 
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PRODUCT OF MATRIX AND INVERSE 

0. 100000E-01 0. 190735E-05 

-0. 707805E-07 0. 999999E-00 

CYCLE NUMBER 1 

NEW CONSTANTS 

UC( 1) = 0. 91269932E-01 

UC( 2) = 0. 97631 855E-01 


CYCLE NUMBER 3 

OLD CONSTANTS 

UC( 1) = 0. 97914543E-01 

UC( 2 ) = 0. 10681258E-02 

INTEGRATION OF PHI AT TIME = 1. 000 
0. 405721E-00 0.419558E-01 

-0. 410815E-00 -0. 42425 9E-01 

SECOND PARTIALS OF S MATRIX 

0.859218E-03 -0. 189624E-03 

-0. 189624E-03 0.445226E-04 

PARTIAL OF S VECTOR 

-0. 1 1 1303E-03 
0. 235987E-04 

INVERSE OF SECOND PARTIALS OF S 

0. 193811E-05 0. 825453E-05 

0. 825454E-05 0. 374026E-06 

DETERMINATE = 0.229721E-08 


MATRIX 
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PRODUCT OF MATRIX INVERSE AND PARTIAL VECTOR 

-0.209212E-00 
-0. 361004E-00 

PRODUCT OF MATRIX AND INVERSE 

0. 999982E-00 -0. 610352E-04 

0. 0 0. 100001E-01 

CYCLE NUMBER 3 

NEW CONSTANTS 

UC( 1)= 0. 10000666E-02 

UC( 2) = 0. 11042262E-02 


CYCLE NUMBER 6 

OLD CONSTANTS 

UC( 1)= 0. 99998312E-01 

UC( 2) = 0. 10999827E-02 

INTEGRATION OF PHI AT TIME = 1. 000 
0.408749E-00 0. 408706E-01 

-0. 408704E-00 -0. 408252E -01 

SECOND PARTIALS OF S MATRIX 

0. 766933E-03 -0. 1 74278E-03 

-0. 174278E-03 0.451110E-04 

PARTIAL OF S VECTOR 

-0. 1 84601E-07 
0. 39741 5E-08 

INVERSE OF SECOND PARTIALS OF S MATRIX 

0. 106786E-05 0.412546E-05 

0.412546E-05 0. 181547E-06 


DETERMINATE = 0. 422444E-08 
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PRODUCT OF MATRIX INVERSE AND PARTIAL VECTOR 

-0. 331756E-04 
-0. 400703E-04 

PRODUCT OF MATRIX AND INVERSE 

0. 999994E-00 0.0 

0. 953674E-06 0. 999998E-00 

CYCLE NUMBER 6 

NEW CONSTANTS 

UC( I)= 0. 99998636E-01 

UC( 2) = 0. 10999867E-02 

CONVERGENCE OF UNKNOWN CONSTANTS 

INTEGRATION OF SYSTEM! USING IDENTIFIED PARAMETERS 
CYCLE NUMiBER 6 
NEW CONSTANTS 


UC( 

1) = 

0. 99998636E-01 

UC( 

2) = 

0. 10999867E-02 

X( 

1, 

1-0) = 

0. 40870443E -01 

X( 

2, 

1-0) = 

-0. 40825028E-01 

X( 

1, 

2-0) = 

0. 55907670E-01 

X( 

2, 

2-0) = 

-0. 55862244E -01 

X( 

1, 

3-0) = 

0. 61439551E-01 

X( 

2, 

3-0) = 

-0. 61394121E-01 

X( 

1, 

4-0) = 

0. 63474596E-01 

X( 

2, 

4-0) = 

-0. 63429177E-01 

X( 

1, 

5-0) = 

0. 64223230E-01 

X( 

2, 

5-0) = 

-0. 64177811E-01 

X( 

1, 

6-0) = 

0. 64498663E-01 

X( 

2, 

6-0) = 

-0. 64453185E-01 

X( 

1, 

7-0) = 

0. 64599991E-01 

X( 

2, 

7-0) = 

-0. 64554513E-01 
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X( 1, 8-0)= 0. 64637244E-01 

X( 2, 8-0)= -0. 64591825E-01 


X( 1, 9-0)= 0. 64650953E-01 

X( 2, 9-0)= -0. 64605534E-01 


X( 1,10-0)= 0. 64655960E-01 

X( 2, 10-0)= -0. 64610541E-01 


4. Hardware and Software Specifications 

(3.) HARDWARE SPECIFICATIONS 

The program was tested on an IBM 360/60 system. The sample 
run presented above uses the values 


(29) 

r = 

2 , 

(30) 

m = 

2 , 

(31) 

N = 

10 , 


and the step used in the integration was h = 0. 01 . The total processor 
time was 28 seconds. Several similar examples were used to test the 
program. 

(b) SOFTWARE SPECIFICATIONS 

The program is written in FORTRAN IV (G). The following table 
lists those parameters that are to be supplied by the user. The param- 
eters are to be punched in cards with.one entry per card consisting of 
an integer i (the variable number)" 1 ' and the value of the variable. The 
table includes the program names of the variable, the symbol used for 


* The variable number has the purpose of helping the user organize his 
data deck by consecutive enumeration of the individual cards on the first 
two columns. It does not serve any other purpose. 
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the variable in the theory (if any), the format in which the card is to be 
punched, and the value of the variable, All the values should be manda- 
torily supplied. Also given in the table is the value of the variable in 
the sample run. 


TABLE 1 

Data Input Requirements 


a) Integer Values 


Program 

Names 

Symbol 

Meaning 

Format 

SampL 

Run 

IC(1) 

- 

run number for use of 
user records 

12,221. 3 

2 

IC (2), NE 

r 

number of linear ordinary 
differential equations 

12, £21. 8 

2 

IC (3 ), NUC 

m 

number of parameters to 
be identified 

12,221. 8 

2 

IC (4), NC3 

k 

number of variables 
observed 

12, 221. 3 

1 

IC (5 ) , NO 

N 

number of observation 
times 

I2,E21. 8 

10 

IC (6), NCYC LE 

-- 

maximum number of 
iterations to be computed 

12, E21. 8 

10 

IC (7), IF LAG 


= 1, observations are to be I2,E21,8 
generated by exact integration 
of x with given values 
“0, this process is omitted 

0 

IC (8), IF LAG1 

" 

= 1, then the matrix 0(t) is 
printed every MP-th step 
= 0, no such printout 

12, E21. S 

1 

IC (9), IF LAG2 


= 1, x(i - 0) at every observa- 
tion point is printed 
= 0, no such output 

12, E21 . 8 

0 
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Program 

Names 

Symbol 

Meaning 

F ormat 

Sample 

Run 

IC (10), IF LAG 3 


= 1, matrix B and vector E 
are printed 
= 0, no such printout 

12, E21. 8 

i 

IC ( 1 1 ) , IF LAG4 

-0 

= 1, B" 1 , det B" 1 , B^E and 
BB - 1 are printed 
= 0, no such printout 

12, E21 , 8 

i 

IC (12), MP 


See IC (8) 


100 

b) Real Variables 




RC (1 ), H 

h 

integration step for integra- 
tion of fundamental 
matrices 

12, E2 1 . 8 

0,01 

RC (2), EPS 

e 

criterion for convergence 

12, E21...8 

0. 001 

RC(3),GAMA(1) 

• • • • • 

Y i 

• • • 

weighting factor for obser- 
vations of first component 
of x 

12, E21 . 8 
• • • 

1.0 

0 9 9 

RC (2+NOS), 
GAMA (NOS) 

\ 

see above 

12, E21. 8 

1.0 

C(l) 

c i 

• » • 

initial condition for 
variable 1 

12, E21. 8 

• • t 

8. 0 

» • 0 

C(NE) 

c 

r 

initial ’condition for 
variable r 

12 , E 2 1 . 8 

8. 0 

UCO(l) 

a i 

* • • 

initial guess for 
parameter 

12, E21. 8 

• • ♦ 

8. 0 

UCO(NUC) 

a 

m 

initial guess for 

parameter a 
m 

12 , E21. 8 

8. 0 
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The option given by IC(7) allows the program to be tested by 
using values produced by exact integration of the given system as experi- 
mental data. This is a good feature to check the reliability and conver- 
gence characteristics of the program for the particular user's system. 


Subroutines Supplied by the Program 


1. INTI and INT2 Adams -Moulton integration routines for the integra- 

tion of the fundamental matrices 0(t), 'f.(t), A. (t) . 

J jq 

2. DAUX Contains the differential equations for the funda- 

mental matrices. Called by INTI and INT2. 

3. MATINV Matrix inversion of B . 


4. MULT ** Storage of matrix B before inversion. Multipli- 

cation of B by B“1 and printout. This routine is 
called if IFLAG4 =1. 


5. OUTCON 

6. OUT1 ** 

7. OUT2 (IB) ** 

8. OUT3 ** 

9. OUT4 ** 


Output of old and new estimates of the unknown 
parameters each cycle. 

Called if IF LAG 1 = 1 to output 0 every MP-th 
step. 

Called if IF LAG 2 = 1 to output the value x(lB-O) 
at time IB . 

Called if IFLAG3 = 1 to output B matrix. 

Called if IF LAG = 1 to output B - ^ , Determ B“^ , 
and B"^E . 


** Fortran printing format is 8E15.6. If more than eight differential 
equations are used, subsequent values are printed in following lines. 
The format may be changed by user. 


Subroutines to be Supplied by User 

1. SYSMAT Input of the system matrix A . Array A (I, J) should 

contain the (I, J)-th element of the system matrix A . 
The equations for the compartments being observed 
must be the first equations of the system. 
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2. PDA 


3. INDOSE 


4. CHANGE 


Input of partial derivatives of A with respect to 
the unknown parameters. Array AP(K,i,J) should 
contain the (I, J)-th element of SA/3a^ . Array 

APP(M, I,J) should contain the (I, J)-th element 

of 5 A/d ada. where M = k (k — 1 )/2 + H for k ^ H . 
k Xj 

Input and printout (if desired) of dosages. DI(I, J) 
should contain dosage to be administered to I-th 
compartment at time J . 

Called if IF LAG = 0 . Input and printout (if desired) 
of observations BI(I, J) should contain the observa- 
tion of the i-th compartment at time J . 


Dimensioning 

To be supplied by user to fit his particular demand and memory 
constraints. 

The following arrays are dimensioned in the main program. 
DIMENSION T(N*10), PD (NUC ) , BB (NUC , NUC ) , IPIV OT (NUC ) , INDEX 
(NUC, 2), PIVOT(NUC ) where N = (1+NUC)*(2+NUC)/2*NE*#2. 

The following COMMON statement with arrays dimensioned as 
indicated should be the first statement of the main program and the sec- 
ond of all subroutines. 

C OMMON DE TERM , EPS , HH, IC YC LE , IF LAG , IF LAG 1 , IF LAG 2 , IF LAG 3 , 
IF LAG4 , JMAX , NN , N2 , NC O N, NC YC LE , NE , NO , NOS , NUC , 

A(NE , NE ), AP(NUC , NE , NE ), APP(M, NE , NE ), B (NUC , NUC ), 
BI(NOS, NO), C(NE), DI(NE , NO),E(NUC , 1 ), GAMA (NOS), 

IC (12), RC(2+NOS), UCO(NUC ), X(NE), XP(NE, NUC ), XPP(NE, M), 

Y (NE ), YP(NE, NUC), YPP(NE, M) 


where M = NUC*(NUC+1 )/2. 



19 


5. Mode of Availability 

Persons interested in obtaining a copy of this program should 

contact 


June Buell 

Dept, of Electrical Engineering 
University of Southern California 
Los Angeles, California 90007 


6. Indexing Information 

System identification, linear systems, drug regimes, pharmaco- 
kinetics. 
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Read Constants 
and Initial 
Approximation 


N = 1 


Integrate Differential Systems to 
Produce Values of the Matrices 

0(1), Y.(l), A. (1) 

J jq 


x(0) = c 

x (0) = 0 

x. (0) = 0 

jq 


E = 0 
n 

B =0 
nq 


Generate x(l), x.(l). 

and x. (1 ) from 

J 

jq 

x(I -1+0), x.(I-l), and x. (I -1) Using 

J 

jq 

Recurrence Relations 

of Step (ii) 


x(l +0) = x(l) +D(1) 


Update Values of E and B by Adding 
I-th Term to Running Sums 


?) 


Solve Linear Algebraic System (Step (iii ) ) 


| a 1 - a 0 I < e 

q q 

For all q? 


Figure 1 

Flowchart of the Method 













